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Abstract 



\ We propose three fast algorithms for solving the inverse problem of the thermoa- 

^ I coustic tomography corresponding to certain acquisition geometries. Two of these 

cH ■ methods are designed to process the measurements done with point-like detectors 

placed on a circle (in 2D) or a sphere (in 3D) surrounding the object of interest, 
g \ The third inversion algorithm works with the data measured by the integrating line 

detectors arranged in a cylindrical assembly rotating around the object. The number 
of operations required by these techniques is equal to 0{n^ logn) and 0{n^ log^ n) for 
j> I the 3D techniques (assuming the reconstruction grid with nodes) and to 0(n^ logn) 

CO ■ for the 2D problem with n x n discretizetion grid. Numerical simulations show that 

our methods are at least two orders of magnitude faster than the existing algorithms, 
without any sacrifice in accuracy or stability. The results of reconstructions from real 
^ I measurements done by the integrating line detectors are also presented, to demonstrate 

O ■ the practicality of our methods. 



Introduction 

Thermoacoustic tomography (TAT) and the closely related photoacoustic modality (PAT) 
are both based on the measurements of acoustic waves generated within the object of in- 
terest by a thermoelastic expansion [17,27]. To initiate the wave, the object is illuminated 
by a short electromagnetic (EM) pulse whose energy is partially absorbed by the tissues. 
The subsequent increase in the local temperature makes the tissues expand, which in turn 
generates the outgoing acoustic wave registered by the detectors surrounding the object. By 
solving the corresponding inverse problem one can reconstruct the distribution of the initial 
pressure inside the body. Spikes in the initial pressure are indicative of cancerous tumors 
that absorb much more EM energy than healthy tissues; thus, cancer detection (particularly, 
in breast imaging) is one of the most promising applications of this modality. General review 
of the inverse problem of TAT/PAT, and investigation of the related mathematical questions 
can be found in reviews [18,35,36] and references therein. 

The recent advances in the measuring technology of TAT/PAT make it possible to collect 
enough data to reconstruct high resolution 3D images of the object of interest. Since such 
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a reconstruction amounts to computing many millions of unknowns, the use of existing 
reconstruction algorithms may lead to inordinate computation time. For example, methods 
based of filtration/backprojection formulas would require several days of computations per 
one high-resolution 3D image. We, thus, propose three fast reconstruction algorithms for 
certain acquisition geometries with spherical or circular symmetry, including, in particular, a 
fast method for the measuring scheme with integrating line detectors arranged in a rotating 
cylindrical assembly. The present algorithms produce high quality images two orders of 
magnitude faster than the existing methods. 

Let us briefly review existing reconstruction techniques while paying close attention to the 
asymptotic estimates of the number of floating point operations (flops) needed to complete 
one reconstruction. (While such estimates hide the unknown constant factor, they usually 
provide a good qualitative measure of the efficiency of the method.) Ideally, one would like 
to reconstruct an image on a grid with unknowns in 0{N) ffops. However, methods that 
require 0(A^log^ A^) ffops (where a is some constant) are still generally considered "fast". 

The simplest inverse problem in TAT/PAT arises when small (point-like) detectors are 
placed on an inffnite plane. In this case the explicit solution can be obtained using the 
Fourier transform techniques (see [3, 10,23] and references therein). The authors of [14], by 
combining proper discretization of such a solution with application of the non-uniform FFT 
(NUFFT) developed a fast algorithm that reconstructs images on n x n x n grid (in 3D) in 
0{n^\ogn) ffops. However, in any practical application such a measuring plane needs to be 
truncated, which leads to a loss of information about the wave fronts nearly parallel to the 
plane. Thus, it is preferable to use closed (and bounded) measuring surfaces. 

One of the simplest closed surfaces is a sphere, and ffrst explicit solutions of the TAT/PAT 
problem in a closed domain were obtained for circular (in 2D) and spherical (in 3D) acquisi- 
tion geometries in [25] and [26] by means of separation of variables. Later, some modiffcations 
of the 2D formulas of [25] were proposed in [2] and [13] in order to avoid the divisions by zero 
in the original formulas of [25]. Such techniques have complexity 0{ii?) ffops for a 2D grid 
with unknowns. A straightforward implementation of the 3D version of the series solu- 
tion [26] results in a rather slow 0{n^) algorithm. However, a modification of this technique 
proposed in [32], if properly discretized, yields a faster 0{n^) algorithm. 

There also exist several explicit inversion formulas of filtration/backprojection type [4, 11, 
12, 19,24,38] for the spherical or circular acquisition geometries. In spite of the theoretical 
importance of these formulas, the corresponding algorithms are not fast, requiring 0{n^) 
fiops for a 2D reconstruction on n x n grid, and 0{n^) flops for a 3D problem. 

The time reversal by means of a flnite difference solution of the wave equation back in 
time [7, 16] is faster (at least in 3D). This technique can be adapted to almost any closed 
acquisition surface, and it allows one to take into account (known) variations in the speed of 
sound within the region of interest (most other techniques are applicable only if the speed 
of sound is a known constant). These methods have complexity 0{n^) in 3D and 0{n^) in 
2D. 

Mathematically equivalent to the time reversal are the methods based on expanding the 
solution of the wave equation into the series of eigenfunctions of the Dirichlet Laplacian in 
the domain surrounded by the acquisition surface [1,20]. This technique is computationally 
efficient only if there exists a fast method for the summation of the eigenfunctions. For 
example, if the detectors are placed on a surface of a cube, such a summation can be per- 
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formed using the 3D Fast Fourier Transform (FFT) algorithm, and one obtains a very fast 
0{n^ logn) reconstruction technique [20]. This is the only known fast method for TAT/PAT 
problems with point-like detectors placed on a closed surface. 

In addition to point-like detectors, there exists another new and promising class of mea- 
suring devices, the so-called integrating hne detectors [5,6,28]. A sensing element in a hne 
detector is a straight optical fiber coupled to a laser interferometer. (In some implemen- 
tations the fiber is replaced by a laser beam propagating through the water in which the 
object of interest is immersed.) Propagating acoustic waves slightly elongate the fiber, and 
this elongation is registered by the interferometer. The measured value is proportional to 
the line integral of the acoustic pressure. Since the fiber can be made very thin, the use 
of such detectors can significantly increase the spatial resolution of TAT/PAT. However, in 
order to enable this new acquisition technique one needs new inversion algorithms. 

The most studied measurement scheme of this sort [5,6,28] utilizes integrating line de- 
tectors arranged in a rotating cylindrical assembly (see Figure 1(a)), with the object placed 
inside the cylinder. A two step inversion procedure was proposed in [28] for solving the 
arising inverse problem (see also [4, 7, 29, 30]). It is based on the observation that, for a fixed 
orientation of the cylinder, the data measured by the line detectors are related to the line in- 
tegrals of the sought function (over the parallel lines) by the two-dimensional wave equation. 
Thus, the latter line integrals can be reconstructed by one of the 2D methods for TAT/PAT 
with conventional point-like detectors placed on a circle. The first step of the reconstruction 
procedure consists in solving such a 2D problem for each orientation of the cylinder. On the 
second step one performs a series of inversions of the 2D Radon transform to reconstruct 
the sought function from the line integrals whose values were found on the first step. Since 
the fastest known algorithms for the circular geometry require at least 0{n^) operations, 
and the number of problems to be solved is (9(n), the first step needs at least O(n^) fiops. 
The second step as described in [4,7,28-30] is also an O(n^) fiops procedure; however, since 
methods for fast ((9(n^logn)) inversion of the 2D Radon transform are well-known (see [22] 
and references therein), it can be accelerated to 0{n^ log n) fiops. Thus, the total number 
of operations (0(n^)) for the whole inversion technique is determined by the first step. 

One of our goals is to develop a fast algorithm for the data acquisition scheme with 
integrating line detectors described in the previous paragraph. To this end, we first develop a 
fast (9(n^ logn) reconstruction algorithm for the 2D problem with point-like detectors located 
on a circle. Since such a technique (and its 3D generalization) are of independent interest 
for problems with conventional detectors, we present them separately in Section 2. Further 
analysis of the 3D problem with the line detectors shows that the two-step reconstruction 
procedure mentioned in the previous paragraph is somewhat suboptimal. Instead, a direct, 
eflScient 0{n^logn) reconstruction algorithm can be built by modifying the 2D method 
described in Section 2.1; we present this technique in Section 3. 

The results of numerical simulations show that our methods yield very fast reconstruc- 
tions without any sacrifice in stability or in the resolution of the images. In addition to 
simulations, in Section 4 we illustrate the work of our algorithms by reconstructing images 
from the data of real measurements performed in RECENDT (Research Center for Non- 
Desctructive Testing, Linz, Austria). 
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(a) 



(b) 



Figure 1: Measuring scheme with hnear integrating detectors (a) general view (b) geometry 

1 Formulation of the problems 

We wiU assume throughout the paper that the speed of sound in the tissue is constant (in 
this case one can set it to be equal 1 without loss of generality), and that the attenuation 
is negligible. This simplified model is acceptable in such important applications as breast 
imaging, and most of the formulas and algorithms mentioned in the Introduction are based 
on these assumptions. (In the situations when the variations in the speed of sound need to 
be taken into account, other techniques (e.g., time reversal) should be used.) 

Under the above assumptions the acoustic pressure u(x^t) solves the following initial 
problem for the wave equation in the whole space 

utt = Au, X e M"^, t G [0,oc), 

^(x,0) = /(x), (1) 
ut{x,0) = 0. 

were the initial pressure of the acoustic wave /(x) is the function we seek, and dimension 
d equals 2 or 3, depending on the context. We will assume that /(x) is finitely supported 
within the region of interest f], and that the measurements of u{x^t) are done outside Vt. 

1.1 Circular geometry with conventional detectors 

As explained in the Introduction, we start with the problem that arises when the pressure 
is measured by the conventional point-like detectors placed on a circle (or, in 3D, a sphere) 
dB of radius i?, and Vt coincides with the corresponding disk (or ball) B. In other words, 
the data P{y^t) are defined as follows 

P{y, t) = u{y, t) \y^Q^ , t G [0, oc). 
Problem 1 (2D) Reconstruct f{x) from P (y, t), x G M^; y G dB, t G [0, oc). 

The solution of this problem will serve as an important step in solving problem 2. 
We will also present a fast algorithm for solving the 3D version of this problem: 
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Problem 1 (3D) Reconstruct f{x) from P {y, t), x eR^, y e dB, t e [0, 2R]. 



The difference in the observation time used in 2D and 3D versions arises since in 3D, due 
to the Huygens principle, u{x^ t) = on dB for all t > 2R. We present fast algorithms for 
solving Problem 1 in 2D and 3D in Section 2. 

1.2 Acquisition scheme with the Une detectors 

Suppose that the acoustic pressure u{x^t) satisfying conditions (1) is measured by the inte- 
grating line detectors lying on the surface of a cylinder of radius R whose axis is parallel to 
the vector D{a) = (cos a, 0, sin a) and passes through the origin (see Figure 1). If we denote 
the left normal to D lying in the horizontal plane (spanned by the vectors ei = (1, 0, 0) and 
63 = (0, 0, 1) ) by N{a) = (— sin a, 0, cos a), then each detector lies on a line /(a, /3) passing 
through the point A{a, /3) = i?cos f3e2 + i?sin f3N{a), where 62 = (0, 1, 0) is the vertical unit 
vector. This detector measures the value proportional to the line integral Pa{y{f3),t) of u : 

P,(y(/5),t) = J u{y,{P)e2 + yMN{a) + sD{a),t)ds, (2) 

y{P) = (i?cos/3,i?sin^). 
where we assume for simplicity that the detectors are infinitely long. 

Problem 2 Reconstruct the initial condition f{x) (supported within a hall of radius R cen- 
tered at the origin) from the measurements Pa{y{l3)^t) known for all a G [0,7r], /3 G [0,27r], 
t G [0,0c). 

We solve this problem in Section 3. 

2 Fast algorithms for Problem 1 

As mentioned in the Introduction, there exist a variety of the methods for the solution of 
Problem 1 in 2D and 3D. However, none of the known methods have optimal computational 
complexity. A fast algorithm for the 3D version of the problem is needed since even relatively 
fast 0{n^) methods require several hours of computing time; methods based on slower 0{n^) 
discrete versions of explicit inversion formulas can easily run a couple of days on larger 
computational grids. 

A single solution of the 2D problem does not take long in practical terms even if a slow 
algorithm is used. However, in order to solve Problem 2, one need to solve the 2D version 
of Problem 1 0{n) times, which can again result in several hours of computation. Hence, a 
fast method is also needed. 
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2.1 2D case 
2.1.1 The algorithm 

Solution u{x,t) to the initial value problem (1) in 2D can be written [34] in the following 
form: 



0, t<\x\. 



u{y,t)^ / f{x)—^{y - x,t)dx 

B 

where ^{x, t) is the free space Green function of the wave equation: 

In particular, the measured data P{y, t) are given by the similar formula 

P{y,t) = J fix)-^^y-x,t)dx, yedB, (3) 

B 

where B is disk of radius R and dB is its boundary. Let us find the Fourier transform P{y^t) 
oiP{y,t) in t : 

P{y, X) = J P{y, t)e''^dt = -i\ J f{x)Hy - A), ye OB, (4) 

M B 

and the Fourier transform A) of the Green function: 

l>(x,A) = J ^{x,t)e''^dt ^ ^H^^\x\x\). 

By combining (3) and (4) we obtain: 

A j f{x)H^'\X\y - x\)dx = 4P(y, A). (5) 

B 

Let us utihze the addition theorem for i/g^^ (see, for example [8]): 

oo 

H^'\x\y -x\)=J2 H\^{XR)J\,\{\r)e^^(^-'\ (6) 

k=—oo 

X = rx, x{9) = (cos ^, sin 
y{R^ (f) = i?(cos (fi^ sin ip) , 
R>r. 

Also, let us expand P{y{R^ (p)^ A) and f{rx{6)) in the Fourier series in p and 6\ 

oo 

k=—oo 
oo 

firm) = E u^y"^'^ (7) 



where coefficients A (A) and /m(^) are given by the formulas 



27r 





27r 



27r/ 





By substituting (6) into (5) and expanding in the Fourier series one obtains 



2tt 

A 



27r oo 



_0 k=-oo 

(27r oo \ 
y J f{rm)M(^r)e-^''rdrdeA (8) 
/ 

and further, by utihzing (7): 

oo 

Pk{X) = ^XHll1{XR) I Mr)J\,\iXr)rdr. 



The above formula relates Pk{X) to the Hankel transform of and, since the latter 

transform is self-invertible. 



oo 



2 f Pk{X) 







Formula (9) was first presented in [2]. A somewhat similar expression was obtained in [25]; 
instead of H^^^^{XR) it contained the Bessel functions J\k\{)^R) in the denominator, and a term 

that corresponds to the real part of our Pfc(A) in the numerator. In that formula, theoretically, 
the zeros of RePfc(A) cancel the zeros of the Bessel function in the denominator. However, 
when Pk{X) is computed from the data of real measurements contaminated by noise, the 
cancellation will not happen automatically. In [13] a technique based on the use of Fourier- 
Bessel series was proposed to avoid such division by zeros. Our approach is based on formula 
(9); since Hankel's functions do not vanish for real (and bounded) values of the argument, 
this formula provides a stable (and simple) way to recover //c(r). 

A straightforward discretization of equation (9) leads to a method that would require 
0{ri^) fiops per each fk and 0{n^) per whole reconstruction, and thus, is not fast enough. 
To develop a fast algorithm we combine (7) and (9) and represent / in the following form: 

/ oo 



firm) = - E / -J^M^'^''"'"^^ I • (10) 
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Let us now consider plain waves W{rx{9)^ A): 

W{rx{9),A) = exp(iAr cos(6^ — ip)), 
A = A(cos(^, sin(^), 

and expand these waves in the Fourier series in Lp. The values of the corresponding Fourier 
coefficients can be found using the Jacobi- Anger expansion [8]: 



so that 



k=—oo 



27T 



By substituting the latter formula in (10) we obtain 

oo 27r r 



f{rx{0)) 







E 



ikif 



.k=—oo \k\ 



W{rx{9),A{X, ip))dipXdX. (11) 



Let us denote by f{A{X^(p)) the expression in the brackets (the choice of such a notation 
will become clear momentarily) 



/(A(A,^)) = - ^ —(1)77^ 



^ikcp 



(12) 



Then (11) can be re- written as 



f{x) = — / /(A)exp(zx. A)dA, 



which means that /(A) is, in fact, the 2D Fourier transform of f{x) defined in the standard 
way: 

/(A) ^ ^ J /(^) exp(-zx • A)dx. 

Formula (12) allows us to compute /(A) for all values of A 7^ 0. In order to find /(O) the 
following integral identity can be used: 



R 



R 



2Po(A) 



00 ^ r R 
f 2Po(A) ? 

J ttH^^\XR) J 



rJo{Xr)dr 



dX 



2Po(A) 
7rXH^^\XR) 



Jo{Xr)dX 



dr 



RJi{XR)dX. (13) 
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Since H^^'' has a logarithmic singularity at A = 0, and Ji has a single root at the latter point, 
the integrand in (13) vanishes as A ^ 0, and the formula is well-defined. Now f{x) can be 
computed by applying the 2D inverse Fourier transform to /(A). We utilize the FFT's on 
various steps of the resulting algorithm to make the computations fast: 
The algorithm for Problem 1 in 2D: 

1. On an equispaced grid in A compute A) using ID FFT in time: 

2. For each value of A in the grid compute Pa;(A) using ID FFT in ip\ 

27r 



A(A) = PiyiR,^),X)e-''^d^. 







3. For each value of A > in the grid compute coefficients 6/c(A) 

,,(A) = ^(-'rp-w, (14) 

4. On a polar grid in A and (p compute /(A(A, (p)) by summing the Fourier series (use ID 
FFT): 

OO 

/(A(A,^))= Yl MA)e^^A^O. 

5. Compute /(O) from formula (13) using the trapezoid rule. 

6. Interpolate /(A) to a Cartesian grid in A. 

7. Reconstruct f{x) by the 2D inverse FFT. 

It is not difficult to estimate the computational costs of the present method. Assuming 
that the computational grid is of size n x n, and the number of detectors and the size of the 
grids in A and ip are 0(n), steps 1, 2, 4, and 7 require (9(n^logn) flops each. Step 3 needs 
0(n^) operations and step 5 is completed in 0{n) flops. 

The interpolation step 6 needs some commentary. It is known [22] that low-order inter- 
polation in spectral domain can be a source of signiflcant error. In particular, interpolation 
in the radial direction (i.e. in A) is more sensitive than that in the angular direction (i.e. in 
(f). The results presented below were obtained by combining linear interpolation in Lp with 
interpolation by cubic polynomials in A. We found that this yields accuracy that is more 
than sufl&cient for practical tomography applications. Since such interpolation is local, the 
associated operation count is 0{n^). However, if a more accurate interpolation is desired, 
one can apply methods based on global trigonometric interpolation using the NUFFT (see, 
e.g. [9]), for the total cost of 0{ri^logn) operations. In any case, the asymptotic estimate 
for the total number of operations required by the present method is (9(n^logn). 
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Figure 2: Reconstruction in 2D (a) phantom (b) reconstruction from noiseless data (c) 
reconstruction in the presence of 50% noise (d) comparison of noiseless (black line) and 
noisy (gray line) data 



2.1.2 Numerical simulations 

In order to evaluate the performance of our algorithm we simulated high resolution projec- 
tions that corresponded to a phantom consisting of several characteristic functions of disks 
lying inside a unit circle, as shown in Figure 2(a). The number of the simulated detectors 
was equal to 272 (this corresponds to the number of detectors in the set of real measure- 
ments presented in Section 4). The detectors were placed on a circle of radius 1.05. The 
measurements were simulated for the time interval [0,5], with time step 0.005 (i.e., 1000 
time samples were simulated). A smooth cut-off was applied at the end of this time interval, 
and the signal for t > 5 was neglected. 

The result of the reconstruction on a grid 1000 x 1000 is shown in Figure 2(b); the 
computation took 0.3 seconds on a desktop computer with a 2.4 GHz Intel Core 2 Duo 
processor. The code was written in Fortran-95, and computations were not parallelized. 

In order to test stability of the present method we added 50% (in norm) noise to 
the simulated projections. Figure 2(c) compares the noiseless and noisy data for the first 
detector. Figure 2(d) presents the result of the reconstruction from the noisy data. One can 
notice low level of noise in the reconstructed image, which is unusual for the inverse problems 
of tomography, where at least mild instability (and hence, some noise amplification) occurs 
almost always. A quick glance at the steps of the algorithm confirms that all the operations 
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(a) (b) (c) 

Figure 3: Comparison with the time reversal using finite differences: (a) image obtained by 
the time reversal; (b) a fragment of the image in part (a); (c) a fragment of the image in 
Figure 2(b) reconstructed by the present method 

are stable: the Fourier transforms and series are isometrics in L^, and in the formula (14) 
Pk{X) is multiplied by a factor that vanishes as A goes to infinity. The reason of such stability 
is the presence of the time derivative in the equation (3) that describes the wave propagation. 
In the absence of this derivative (i.e. for example if the initial value of the pressure was zero 
but the time derivative was non-zero) the stability properties of the problem would be similar 
to those of the classical inverse Radon transform in 2D, and some noise amplification would 
occur during the reconstruction. 

It is also interesting to compare the present technique with the results obtained by the 
time reversal computed using finite differences. In order to do so, we solved the 2D wave 
equation back in time in the unit square Jl, using the explicit leapfrog scheme 

(At)2 (Ax)2 

starting with u — 0^ -^u = at time t = 5, and enforcing the boundary condition u[y^t) — 
P{y,t) for all y lying on the boundary dVt of Vt. We chose to do the reconstruction in a 
square domain since the exact boundary conditions are easy to enforce in the nodes of 
the computational grids lying on the sides of the square. (For a domain of a different 
shape enforcing the boundary conditions requires interpolation (see e.g. [7]) which brings 
additional error). In order to guarantee the stability of this explicit scheme we had to increase 
the number of time steps from 1000 to 3800. 

The result of the reconstruction by the time reversal is shown in Figure 3(a). Larger 
features of the image are reconstructed quite well. However, the second-order finite difference 
scheme we used is not very accurate on higher spatial frequencies, which leads to a typical 
diamond-shaped distortion of small round shapes, clearly visible in the magnified fragment 
of the image presented in Figure 3(b). The same part of the image in Figure 2(b) (obtained 
using the present algorithm) is shown in Figure 3(c). The artifacts are much smaller, in spite 
of the lesser number of time steps in the data. 

The time reversal took 214 seconds on the same computer, as before. Even if one discounts 
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this number by a factor of 4 (to account for a larger number of the time steps), our technique 
is still faster by two orders of magnitude. 

2.2 3D case 

In this section we briefly outline the 3D version of the method presented in Section 2.1.1. In 
3D, the Green's function for the wave equation satisfying radiation condition at infinity has 
the following form [34]: 

, 5{t-\x\) 
^ ^ 4:7r\x\ 

As in 2D case, we compute the Fourier transform in t from the data t) 

P{y, X) = J P{y, t)e''^dt = -iX J f{x)Hy - x. A), ye OB, (15) 

M B 

and the Fourier transform $(x, A) from the Green's function 

By combining (15) and (16) we obtain: 

A' J f{x)h^o\My - ^\)dx = 47rP(?/, A). (17) 

We will utilize spherical harmonics Yj^{z)^ z G S^, normalized so that 



where 5k,p is the Kronecker symbol. Let us extend f{x) and P(y, A) in spherical harmonics: 



f(rx) ^J2Y1 fsAr)yi'i^),^^ S^r = |x|, (18) 

5=0 p= — S 

fsAr) = J f{xT)YF^dx, 

oo s 

P(i?y,A) = ^^ A,p(A)17(y), (19) 

5=0 p= — S 

A,p(A) = J P{Ry,X)Yfiv)dy. (20) 

S2 
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As in the 2D case, we make use of the addition theorem for /iq^^ [8]: 

oo k 

h^\\\z - x\) = At^Y. E ht\\\z\)3,{\\x\)Y^{z)Y^l 

k=0 m=—k 

where jk{-) and h^j^\-) are, respectively, the spherical Bessel and Hankel functions. 
By substituting the above equation, together with (18), into (17) we obtain 



oo / 
« /> / oo s 

^(?/,A) = J J (X^5^/.,p(r)y/(x) 



§2 

OO k 



s=0 p=—s 



oo k 



E hf\XRMXr)Y,-{y)Y,-{x) 



.k=0 m=—k 



dxr dr 



= E E Uk\^m'J fkAr)jk{Xrydr\Y, 



k \y)^ 



k=0 m=—k 



and further, by comparing with (19): 

4p(A) = A2/iW(Ai?) 



(OO 
//., 




,p{r)js{Xr)r dr . 



(21) 



The spherical Bessel functions are related to their cylindrical counterparts by the equation 



js{t) = \I^Js^l/2{t). 

Thus, Ps,p{X) can be expressed in terms of the Hankel transform of y/rfs^p{r): 

4p(A) = ^X'/'hi'\XR) (j[V^fsAr)] Js^i/2{Xr)rdr^ . 
Since the Hankel transforms are self-invertible, one recovers fs,p{j^) as follows: 

oo ^ oo ^ 



This approach is close (although not quite identical) to the solution obtained in [26]. 
(An expression equivalent to (22) was derived in [32].) A straightforward computation of 
(22) for all s < \p\ < s leads to an algorithm of complexity at least 0{n^) for an n x 
n X n computational grid. In order to accelerate the computations we choose the following 
approach. By substituting (22) into (18) we obtain: 



{ s=0 p=-s i^s [XK) 



(23) 



s=0 p=—s ^ 
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A convenient integral representation for the term js{Xr)Y^(x) in the above equation is given 
by the Funk-Hecke formula [8]: 



Combining (23) and (24) yields: 



fix) 




EE 



(27r)3/2 / 



EE V:;; 



5=0 p= — S 



A2/i(^\|A|i?) ' V|A| 



where A = zX. Let us denote the term in the brackets by -F(A): 

oo s nr 

Then 



s=0 p=—s 



2z«A,p(i?,|A|)^,/AV 



A^h^^\\A\R) ' V|A|y 



(24) 



(25) 



(26) 



and it becomes clear that function F{A) is the 3D Fourier transform of f{x) defined as 
follows: 



Formulas (15) and (20) allow us to reconstruct the Fourier transform F{A) for all A 7^ 0. 
In order to find F(0) we notice that 



F(0) = 



(27r)3/2 

where /o,o is given by (22), so that: 



f{x)dx 



1 



R 



(27r)3/2 



m = 



1 



V2(7r)5/2 7 h^^xR) 



jo{Xr)dXr'^dr 



Po,o{R, A) 
V2(7r)V2 7 h^^\XR) 



R 



jo{Xr)r'^dr 



dX. 
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Since jo{t) — sin(t)/t, the integral in the brackets can be easily evaluated: 

R 

A2 V \R 





leading to the following expression for F(0): 



J joiXrydr = ^ f ^^^^ - cos Ai? ) , 



OO ^ 

R f Poo(RA) /sinAi? , \ 

F(0) = / "'"\ ' ' — cos XR dX 

OO 

-Po,o{R, A) exp(-iAi?) f cos XRj d\. 



V2(7rf/^J A 



Thus, f{x) can be reconstructed by the following method: 
The algorithm for Problem 1 in 3D: 

1. On an equispaced grid in A compute P(y, A) using ID FFT in time: 

P(y,A) = J P{y,t)e^'^dt 

2. For each value of A in the grid expand P{Ry^ A) in spherical harmonics in 

3. For each value of A > in the grid compute coefficients bs^p{X): 

f2 i-P.,,(R,\) 

4. On a spherical grid in A compute F(A) by summing spherical harmonics: 



OO s / A \ 

s=0 p=-s ^ ' ' ^ 



5. Compute F(0) from formula (27) using the trapezoid rule. 

6. Interpolate F{A) to a Cartesian grid in A. 

7. Reconstruct f{x) by the 3D inverse FFT. 
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1 1 

Figure 4: 3D phantom 

The most time-consuming steps of this algorithm are steps 2 and 4, corresponding to 
the Fourier analysis and synthesis on a sphere. One of the simplest ways to decompose a 
function in the spherical harmonics combines expansion into the Fourier series (which can 
be done fast using the FFT) with the expansion into the series of Legendre functions. The 
straightforward implementation of the latter Legendre transform requires 0{ri^) flops for a 
one dimensional function deflned by n samples. Then, expanding one 2D function on a sphere 
(deflned by n x n samples ) into a series of spherical harmonics needs 0{n^) operations, and 
the present method involving 0{n) of such expansions would require O(n^) flops (synthesis 
of spherical harmonics is quite similar to the analysis). 

Recently, several versions of the Fast Spherical Harmonics (FSH) transform were intro- 
duced by several research groups [15,21,31,33]. FSH is asymptotically fast; it performs 
Fourier analysis and synthesis on a sphere in (9(n^log^n) flops. If this algorithm is used on 
steps 2 and 4 of our reconstruction technique, the resulting method will also be fast, requiring 
0{n^ log^ n) operations per a 3D reconstruction. One has to keep in mind, however, that the 
break-even size for which the FSH starts to outperform the simple slow method described 
in the previous paragraph is currently of order of hundreds, and thus our 3D reconstruction 
method would not signiflcantly outperform slower (9(n^) techniques for the sizes of the prob- 
lems we consider in this paper. In addition, ready-to-use implementations of the FSH are 
only available for a restricted number of operating systems and computational languages. 
However, as the work on the FSH progresses, the performance and the ease of programming 
of our 3D method will improve. 

3 Fast algorithm for Problem 2 

It is known [4,28] that Problem 2 can be reduced to solving a set of Problems 1 in 2D, 
followed by a set of numerical inversions of the 2D Radon transform. Indeed, let us fix unit 
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vector D(a) and consider all line integrals Va{h^t) in the direction D(a) of u{x^t) : 

Va{h^ t) — J u{hiN{a) + /i2e2 + sD{a)^t)ds^ 

Ml 

h = (/li,/l2). 




(a) Plane X3 = —0.5 



(b) Plane X2 = -0.5 



(c) Plane xi = —0.5 



Figure 5: Simulation in 3D; first row represents the phantom; second row is the recon- 
struction from noiseless data; third row shows reconstruction from the data with added 50% 
noise 

It is well-known (see e.g. [34]) that integrals Vo^{h, t) satisfy the wave equation in (the 
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"method of descent" is based on this fact): 

Q2 Q2 Q2 



dhl dhl df" 

with the initial conditions 

Va{h, 0) = Ma{h) = J u{hiN{a) + /i2e2 + sD{a), 0)ds 

= J f{hiN{a) + /i2e2 + sD{a))ds, 

In other words, the initial values Mo, of are the integrals of f{x) along the straight 
lines parallel to the vector D{a). On the other hand, integrals Pa{y{f3),t) measured by the 
integrating line detectors (see equation (2) are equal to the values of Va{h^t) for all h lying 
on the circle dB of radius R centered at the origin /i = 0: 

Va{y, t) = P^(i/, t), y^dB, te [0, oc). 

Therefore, for each fixed angle a, values of Ma{h) can be recovered within the disk B 
(bounded by dB) by applying the 2D algorithm of Section 2.1 to Pa{y^t). Further, for a 
fixed value of the vertical coordinate /i2, the values of Mq,((/ii, /12)) correspond to the 2D 
Radon transform of the restriction of /(x) to the plane X2 — Thus f{x)\^^^^^ can be re- 
constructed from Mc^((/ii, /12)) using one of the well-known inversion algorithms for the latter 
transform (see e.g. [22]). In particular, in [4,28] the well-known filtration/backprojection 
algorithm is utilized for such a computation. 

The above mentioned techniques are not as fast as we would like, however. In order to 
accelerate the computation, one can combine the fast method proposed in Section 2.1, with 
the fast Fourier transform-based inversion of the Radon transform [22]. Notice, however, 
that such two stage approach adds to the algorithm a numerical inversion of the 2D Radon 
transform, which may increase the computational error, since such inversion is a (mildly) 
ill-posed problem. Moreover, such an algorithm is somewhat redundant. Indeed, on the 
first stage the values of Ma{h) are reconstructed by means of the Fourier synthesis (step 7 
of the algorithm presented in Section 2.1.1). On the second stage, as a part of the Fourier 
reconstruction from projections, the ID Fourier transform in hi of projections Mc((/ii, /i2)) 
is computed for each value of a 

This redundancy can be eliminated as follows. For each fixed a let us apply the first 5 
steps of the algorithm from Section 2.1.1 to the data Pa{y^ t) and compute functions /q;,2d(A), 
where the subscript indicates the dependence of / on the angle a, and the fact that /q;,2d(A) 
are solutions of 2D problems. For a fixed a, function /a,2D(A) is the 2D Fourier transform 
of Ma{h). Since Ma{h) are the X-ray projections of the initial condition /(x), x G in the 
direction D(a)^ by the well-known slice-projection theorem [22] /a,2D(A) coincides with the 
values of the 3D Fourier transform fsoi^) of f{x) restricted to the plane normal to D{a) 
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and passing through the origin A = 0. Therefore, fsoi^) can be reconstructed directly from 
the values of fa,2D{^) computed for all a from to tt. The sought function f{x) is then 
obtained by computing the 3D inverse Fourier of f^jj^A). This can be summarized in the 
form of the following 

Algorithm for solving Problem 2 

1. For all a from to tt compute fa,2D{^) from Pa{y^ t) by performing the first 5 steps of 
the algorithm from Section 2.1.1 

2. Interpolate values /a,2D(A) to obtain /3d(A) on a 3D Cartesian grid 

3. Compute f{x) from fsoi^) by using the 3D inverse FFT. 

If the number of the cylinder directions is of order of step 1 of the above algorithm 

requires 0{n^logn) flops. The second step involves interpolation in the spectral domain 
from spherical grid to the Cartesian grid. In our simulations we utilized cubic interpolation 
in the radial direction and bilinear interpolation in the angular directions; the number of 
operations involved in this step is 0{n^). The third, flnal step requires 0{n^ logn) operations; 
the whole algorithm then needs only 0{n^logn) flops. 

In order to evaluate the performance of the present algorithm we simulated the measure- 
ments corresponding to 272 integrating line detectors rotated over 512 directions a equis- 
paced from to tt. As a phantom we used a collection of characteristic functions of balls 
centered on the (pair-wise) intersections of the planes Xi — —0.5, X2 = —0.5, Xs — —0.5, and 
lying within the unit sphere, as shown in Figure 4. (The orientation of the axes corresponds 
to vectors ei, 62 and 63 as illustrated in Figure 1(b).) The cross-sections of the phantom are 
shown in the flrst row of Figure 5. 

For each direction a, 500 time samples were simulated in the interval t G [0, 5]; the rest 
of the signal was neglected. The image was reconstructed on a 500 x 500 x 500 Cartesian 
grid containing 125 million unknowns (although only about a half of them lied within the 
unit sphere where function f{x) was supported). The cross sections of the reconstructed 
/(x) are shown in the second row of Figure 5. 

The above computation took 67 seconds on the desktop computer described in Section 
2.1.1. A comparison can be made with a reconstruction obtained by the time reversal 
using finite diflFerences in a cubic domain with 251 x 251 x 251 computational grid. On 
our computer it took about 50 min. Since such a time reversal method scales as (9(n^), 
on a grid of the size 500 x 500 x 500 the reconstruction would take about 13 hours, or 
three orders of magnitude longer than the time required by the present algorithm. (Such 
a comparison is quite crude since a different problem is solved in a different computational 
domain; nevertheless it indicates that our method is indeed very fast.) 

The third row of Figure 5 demonstrates images reconstructed from the data with added 
simulated noise with intensity 50% (in norm) of the signal. The level of noise in the 
reconstructed images is surprisingly low. The explanation of such stability is the same as in 
Section 2.1.1. 
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(a) 



(b) 



Figure 6: Reconstruction of a 2D projection from the data measured by integrating line 
detectors (a) the data (b) the reconstruction 

4 Applications to real data 

In this section we illustrate the work of our algorithms by processing a set of real data kindly 
provided to us by RECENDT. This set of data was measured by integrating line detectors, 
as shown in Figure 1 and described in detail in Section 3. As the test object the researchers 
from RECENDT used a piece of a human hair tied in a knot. The number of the detector 
directions D{a) in this set was 25. For each direction, the linear detector was respectively 
placed in each of 272 equispaced positions on the surface of the (imaginary) cylinder. (In 
fact, 11 out of 272 positions at the bottom were unavailable and the corresponding data 
were replaced by zeros.) For each position 10000 time samples were measured. To reduce 
the noise, the signal was smoothed by a convolution with a Gaussian and downsampled by a 
factor of 10, so that the number of time samples on the input of the reconstruction algorithm 
was 1000. Moreover, to reduce strong artifacts at the beginning and at the end of each time 
series, the signal was set to zero there as well. The resulting set of data (corresponding to the 
detectors being aligned along 63) is represented by the gray-scale image shown in Figure 6(a). 
Each vertical line in this image corresponds to the time series for one detector position of a 
line detector; the bottom corresponds to t = 0. 

Figure 6(b) demonstrates the image reconstructed by applying algorithm of Section 2.1.1 
to the data shown in Figure 6(a). As explained in Section 3, the result is not the complete 
reconstruction, but the X-ray transform of the 3D density function f{x) corresponding to 
the hair. 

In order to reconstruct /(x) we applied algorithm of Section 3 to the full set of data. The 
results are shown in Figure 7. The three images in this figure correspond to the horizontal 
cross-sections of the test object at the levels A, B, and C, respectively. One can see sharp 
spikes in the images corresponding to the location(s) of the hair. There are significant radial 
artifacts in these images, arising due to the insufficient sampling in angle a (i.e. insufficient 
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Figure 7: 3D reconstruction from real data; slices parallel to the plane OxiXs at the X2 levels 
that correspond, respectively, to marks A,B, and C in Figure 6 

number of the detector directions). It is known [22] that for the inversion of the 2D Radon 
transform the optimal number of the angular directions is of the same order of magnitude as 
the desired resolution of the spatial grid (several hundred, in out case). Unfortunately, our 
data set only contained data corresponding to 25 detector directions. As the reconstructions 
obtained in Section 3 from simulated data show, if the number of directions equals several 
hundred, the algorithm yields very accurate and detailed images. 

Acknowledgements 

The author would like to thank RECENDT for providing the data of real measurements. We 
are grateful to Drs. P. Burgholzer, H. Griin, and H. Roitner for the fruitful discussions of the 
data acquisition schemes with integrating line detectors, related data processing techniques 
and open problems in this area. The author gratefully acknowledges support by the NSF 
through the grant DMS-090824. 

References 

[1] M. Agranovsky and P. Kuchment, Uniqueness of reconstruction and an inversion proce- 
dure for thermoacoustic and photoacoustic tomography with variable sound speedy Inverse 
Problems, 23:2089-102, 2007. 

[2] G. Ambartsoumian and P. Kuchment, A range description for the planar circular Radon 
transform SIAM J. Math. Anal, 38(2): 681-92, 2006. 

[3] L. E. Andersson, On the determination of a function from spherical averages SIAM J. 
Math. Anal., 19(1): 214-32, 1988. 

[4] P. Burgholzer, J. Bauer-Marschallinger, H. Griin, M. Haltmeier, and G. Paltauf, Tem- 
poral back-projection algorithms for photoacoustic tomography with integrating line de- 
tectors^ Inverse Problems, 23:S65-S80, 2007. 



21 



[5] P. Burgholzer, C. Hofer, G. Paltauf, M. Haltmeier, O. Scherzer, Thermoacoustic to- 
mography with integrating area and line detectors^ IEEE Transactions on Ultrasonics, 
Ferroelectrics, and Frequency Control, 52(9):1577-83, 2005. 

[6] P. Burgholzer, C. Hofer, G. J. Matt, G. Paltauf, M. Haltmeier, and O. Scherzer, Ther- 
moacoustic tomography using a fiber-based Fabry- Perot interferometer as an integrating 
line detector, Proc. SPIE 6086, 434 42, 2006. 

[7] P. Burgholzer, G. J. Matt, M. Haltmeier, and G. Paltauf, Exact and approximative imag- 
ing methods for photoacoustic tomography using an arbitrary detection surface, Phys 
Review E, 75, 046706, 2007. 

[8] D. Colton and R. Kress, Inverse acoustic and electromagnetic scattering theory, 
Springer- Verlag, 1992 

[9] A. Dutt and V. Rokhlin, Fast Fourier Transforms For Nonequispaced Data, Siam J. Sci. 
Comput., 14(6): 1368-93, 1993. 

[10] J. A. Fawcett, Inversion of n- dimensional spherical averages, SIAM J. Appl. Math., 
45(2): 336-41, 1985. 

[11] D. Finch, M. Haltmeier and Rakesh, Inversion of spherical means and the wave equation 
in even dimensions, SIAM J. Appl. Math., 68(2): 392-412, 2007. 

[12] D. Finch, S. Patch and Rakesh, Determining a function from its mean values over a 
family of spheres, SIAM J. Math. Anal., 35(5): 1213-40, 2004. 

[13] M. Haltmeier, O. Scherzer, P. Burgholzer, R. Nuster, and G. Paltauf, Thermoacoustic 
Tomography And The Circular Radon Transform: Exact Inversion Formula, Mathemat- 
ical Models and Methods in Applied Sciences, 17(4): 635-55, 2007. 

[14] M. Haltmeier, O. Scherzer and G. Zangerl, A Reconstruction Algorithm for Photoacous- 
tic Imaging Based on the Nonuniform EFT, IEEE Trans. Med. Imag., 28(ll):1727-35, 
2009. 

[15] D. M. Healy Jr., D. N. Rockmore, P. J. Kostelec, and S. Moore. EFTs for the 2-Sphere 
- Improvements and Variations, J. Fourier Anal, and Appl., 9(4): 341-85, 2003. 

[16] Y. Hristova, P. Kuchment, and L. Nguyen, On reconstruction and time reversal in ther- 
moacoustic tomography in homogeneous and non-homogeneous acoustic media, Inverse 
Problems, 24: 055006, 2008. 

[17] R. A. Kruger, P. Liu, Y. R. Fang, and C. R. Appledorn, Photoacoustic ultrasound 
(PAUS) reconstruction tomography, Med. Phys., 22: 1605-09, 1995. 

[18] P. Kuchment and L. Kunyansky, Mathematics of Photoacoustic and Thermoacoustic 
Tomography, in Handbook of Mathematical Methods in Imaging, Scherzer, Otmar (Ed.) 
Springer, 2011. 



22 



[19] L. Kunyansky, Explicit inversion formulae for the spherical mean Radon transform^ 
Inverse Problems, 23: 737-783, 2007. 

[20] L. Kunyansky, A series solution and a fast algorithm for the inversion of the spherical 
mean Radon transform^ Inverse Problems, 23: S11-S20, 2007. 

[21] M. J. Mohlenkamp, A Fast Transform for Spherical Harmonics^ J. Fourier Anal. Appl., 
2: 159-84, 1999. 

[22] F. Natterer, The mathematics of computerized tomography, New York, Wiley, 1986. 

[23] F. Natterer and F. Wiibbeling, Mathematical Methods in Image Reconstruction, ser. 
Monographs Math. Model. Comput., Philadelphia, PA: SIAM, 2001, vol. 5. 

[24] L. Nguyen, A family of inversion formulas in thermoacoustic tomography^ Inverse Prob- 
lems and Imaging, 3(4): 649-675, 2009. 

[25] S. J. Norton, Reconstruction of a two-dimensional reflecting medium over a circular 
domain: exact solution^ J. Acoust. Soc. Am., 67: 1266-1273, 1980. 

[26] S. J. Norton and M. Linzer, Ultrasonic reflectivity imaging in three dimensions: ex- 
act inverse scattering solutions for plane, cylindrical, and spherical apertures^ IEEE 
Transactions on Biomedical Engineering, 28: 200-202, 1981. 

[27] A. A. Oraevsky, S. L. Jacques, R. O. Esenaliev, and F. K. Tittel, Laser-based ptoacoustic 
imaging in biological tissues, Proc. SPIE, 2134A:122-128, 1994. 

[28] G. Paltauf, R. Nuster, M. Haltmeier, and P. Burgholzer, Thermoacoustic Computed 
Tomography using a Mach-Zehnder interferometer as acoustic line detector, Appl. Opt., 
46(16):3352-8, 2007. 

[29] G. Paltauf, R. Nuster, M. Haltmeier and P. Burgholzer, Experimental evaluation of 
reconstruction algorithms for limited view photoacoustic tomography with line detectors, 
Inverse Problems, 23: S81-S94, 2007. 

[30] G. Paltauf, R. Nuster, and P. Burgholzer, Weight factors for limited angle photoacoustic 
tomography, Phys. Med. Biol., 54: 3303-14, 2009. 

[31] D. Potts, G. Steidl and M. Tasche, Fast and stable algorithms for discrete spherical 
Fourier transforms.. Linear Algebra Appl., 275/276: 433-450, 1998. 

[32] A. G. Ramm, Injectivity of the spherical means operator, C. R. Math. Acad. Sci. Paris, 
335(12): 1033-38, 2002. 

[33] R. Suda and M. Takami, A fast spherical harmonics transform algorithm. Mathematics 
of computation, 71(238): 703-15, 2001. 

[34] V. S. Vladimirov, Equations of mathematical physics. (Translated from the Russian by 
Audrey Littlewood. Edited by Alan Jeffrey.) Pure and Applied Mathematics, 3 Marcel 
Dekker, New York, 1971. 



23 



[35] L. Wang, (Editor), Photoacoustic imaging and spectroscopy, CRC Press, Boca Raton, 
FL, 2009. 

[36] L. V. Wang and H. Wu, Biomedical Optics. Principles and Imaging, Wiley-Interscience, 
2007. 

[37] M. Xu and L.-H. V. Wang, Time-domain reconstruction for thermoacoustic tomography 
in a spherical geometry, IEEE Trans. Med. Imag., 21: 814-822, 2002. 

[38] M. Xu and L.-H. V. Wang, Universal back-projection algorithm for photoacoustic com- 
puted tomography, Phys. Rev. E, 71:016706, 2005. 



24 



